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Abstract 

Given a next-to-leading order calculation, we show how to set up a computer 
program, which generates a sequence of unweighted momentum configurations, each 
configuration containing either n or n+1 four-vectors, such that for any infrared safe 
observable the average over these configurations coincides with the NLO calculation 
up to errors of order i/res- The core of the algorithm is a method to combine real 
emission and virtual corrections on a point-by-point basis in hard phase space. The 
algorithm can be implemented on top of existing NLO calculations. 



1 Introduction 



As high energy experiments move towards a new generation of colhders, accurate sim- 
ulation of coUision events becomes increasingly important. Among the main tools used 
today are next-to-leading order (NLO) calculations and event generators. Up to now 
these two tools have a complementary regions of applicability. A NLO-calculation for 
an n-jet observable calculates this observable as a sum of two terms, one term living 
on the n-parton phase space, the other on the (n + l)-parton phase space. These two 
contributions correspond roughly to the virtual corrections and the real emissions, where 
infrared divergences have been taken care of by a method like phase space slicing or the 
subtraction method. Usually the two contributions to an observable are only combined 
in the end. As a side-effect, one deals at intermediate steps with variable and/or negative 
weights. 

Event generators like HERWIG |I| and PYTHIA on the other hand produce unweighted 
events through a three-stage process: From an initial hard process parton showers are gen- 
erated. A hadronization model converts at a low scale the partons into hadrons. The hard 
process is usually based on Born matrix elements. 

The two approaches have complementary strenght: Parton showers describe quite rea- 
sonably multiple coUinear emission, whereas a NLO-calculation is by far more accurate 
when n partons are well separated in phase space. 



It is a topic of current interest to merge the benefits of the two approaches into an "event 
generator accurate at NLO". Friberg and Sjostrand have suggested a possible road in 0. 
Similar approaches have been considered by Collins IQ, Potter [|5[ and Dobbs and Lefeb- 
vre and applied to deep- inelastic scattering (0, @) and the production of electroweak 
gauge bosons in hadron-hadron collisions (|^, [[10]). Mrenna has considered the case how 
to combine resummed calculations with parton showers Despite this progress, there 
is still room for improvement: The algorithm of Collins may generate events with nega- 
tive weight. The algorithm of Potter uses phase space slicing and introduces a unresolved 
region, such that the sum of the Born term, the virtual corrections and the real emission 



part is exactly zero |[T^. It is not guaranteed that this region is always small enough such 
that the approximations of phase space slicing are still valid. 



In combining parton showers with NLO-calculations two problems can immediately be 
identified: Negative weights and double counting. In this paper we would like to address 
the first problem. As already stated before, a specific momentum configuration within a 
NLO calculation might come with a negative weight. There is nothing wrong with this 
fact by itself and one might think of carrying this weight along the further simulation of 
the event. However, there are some concerns raised against this approach: 

a. ) Some hadronization models have problems with negative weights 0. 

b. ) Experimentalists are unwilling to accept negative weights. Event reconstruction and 

detector simulation are quite expensive in terms of computer time and an event with 
a negative weight would "annihilate" parts of the CPU-time already spent. 
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Therefore it is desirable to avoid negative weights from the very beginning. In view of 
point b.) this is more a question of efficiency than of principle. 



In this paper we consider an "NLO partonic event generator". Given a next-to-leading 
order calculation, we show how to set up a computer program, which generates a sequence 
of unweighted momentum configurations, each configuration containing either n or n + 1 
four-vectors, such that for any infrared safe observable the average over these configura- 
tions coincides with the NLO calculation up to errors of a resolution variable Hres- The 
basic idea is simple: Consider a (n + l)-parton configuration. If all partons are separated 
by at least yres, we return this (?t, + l)-parton configuration. If on the other hand a parton 
is not resolved at a scale i/res, we average over the real emissions up to the scale i/res, com- 
bine the result with the virtual corrections and return a n-parton configuration. In this 
paper we show that this can be done in a general and effective way: Our procedure does 
not depend on the specific hard process and the integration over the real emission part 
involves only a three-dimensional integration. These two points are important: Today 
many NLO-calculations already exist and we would like to be able to use them without 
major modifications. Our algorithm can take as an input subroutines, which correspond 
to the matrix elements for real emission and virtual corrections without modifications. 
Only the subtraction terms within the dipole formalism have to be implemented. Sec- 
ondly, if a parton is not resolved we have to integrate out this contribution. An important 
ingredient is a factorization of the phase space into a "hard" n-parton configuration and 
the insertion of a "soft" particle into this "hard" configuration. This allows us to inte- 
grate out the unresolved region on a point-by-point bases in hard phase space. We discuss 
explicitly the case of massless QCD processes in electron-positron annihilation, although 
the extension to initial-state partons and/or massive partons is in principle feasible. 

In this paper we do not consider how to attach parton showers to the events nor do 
we address the problem of double counting. This is a separate problem and addressed for 
example in [Q] and [Q. 

This paper is organized as follows: The next section gives a more detailed outline of 
the algorithm. Section 3 reviews the structure of infrared divergences and the dipole for- 
malism. Section 4 deals with basic Monte Carlo techniques and the generation of phase 
space. In section 5 we present an algorithm to generate unweighted events according to 
leading-order matrix elements. In section 6 this algorithm is generalized to unweighted 
events according to next-to-leading order matrix elements. Finally, section 7 contains the 
conclusions. 

2 Preliminaries 

We assume that we have at our disposal a NLO-calculation, consisting of the virtual 
corrections, symbolically denoted by 



da 



V 




) + {M 




1-^ Born) 5 




the Born amplitude for the emission of an additional parton 




(2) 
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as well as colour-ordered partial amplitudes for the Born process. The partial amplitudes 
are needed for the dipole formalism. We denote the Born matrix element by 

dcT^ = {M^sLlMP^rn)- (3) 

We distinguish three scales: 

• The scale ycut- This scale enters directly the definition of a physical observable. An 

example would be the measurement of the 3-jet cross section in electron-positron an- 
nihilation with the Durham algorithm at ycut = 0.01. Results of physical observables 
may depend on the scale ycut- 

• The scale yres'- The resolution scale i/^es gives the lower bound to which we might 
reasonably push ycut- On the experimental side yres is set by the finite resolution of 
the detector. We are here concerned with fixed-order perturbation theory. Consider 
first the situation where we would like to generate a sequence of n four-vectors 
according to the leading-order matrix element. This matrix element is divergent 
when one parton becomes either soft or collinear. We would therefore generate 
predominately events where one parton is soft or collinear. For any infrared-safe 
n-jet observable these events will be rejected since they do not pass the jet-defining 
cuts. This is inefficient. We therefore introduce a resolution function 0„ and 

yres 

generate only events which are at least resolved at the scale yres- Let us now 
consider next-to-leading order calculations. With multiple scales in the problem, 
this will inevitably lead to logarithm of ratios of these scales. If the various scales 
are not of the same magnitude, these logarithm can become large and spoil the 
validity of a fixed-order perturbative expansion. Consider the following situation 
in the real emission part of an NLO-calculation: Assume that we have n partons 
well separated in phase space with an additional parton close to one of the first n 
partons. We would like to replace this (n + l)-parton configuration with an n-parton 
configuration, where the unresolved parton has been integrated out. Integrating out 
the unresolved parton in the real emission part up to a resolution scale yres will yield 
a positive, but divergent contribution. The divergent part is compensated by the 
virtual corrections, and the sum of the real emission part and the virtual corrections 
is finite, but not necesarrily positive. The sign will depend on the scale yres- If yres 
is chosen too small, we do not include enough contributions from the real emission 
amphtude and the combined result will be negative. This is unacceptable for a 
probabilistic interpretation. We therefore have too chose y^es large enough that for 
a given "hard" n-parton configuration we can integrate out an additional "soft" 
parton, which is within yres of one of the "hard" partons and the combination of the 
obtained real emission part with the virtual corrections gives a positive contribution. 

• The scale ymin'- This is an internal technical parameter we introduce for integrating 
out efficiently unresolved contributions. In some place we will neglect contributions 
of order ymin, but since ymin « Ures the introduced systematic error is negligible 
against the approximations made in the previous step. 

Of course the technical problem is to integrate out efficiently the unresolved region for a 
given "hard" n-parton configuration. We show that this involves only a three-dimensional 
integration. We take special care that this can be done efficiently: Our algorithm is based 
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on the subtraction method. This is more efficient than the phase space shcing approach, 
since we don't have to reproduce numerically a logarithm of Umin) which cancels in the 
end. Furthermore we use a specific remapping of the unresolved phase space, which tends 
to fiatten the integrand. 

Our method relies on the following ingredients: 

• The structure of the infrared divergences is universal. We may therefore add and 
subtract an appropriately chosen term, such that the integral over the real emission 
part in the unresolved region and the virtual corrections become finite. Explicitly 
we are using the dipole formalism of Catani and Seymour [0 together with slight 
modifications given by Nagy and Trocsanyi . 

• An algorithm to generate the "hard" n-parton phase space from a set of random 
numbers distributed in [0,1]. There are several algorithms on the market and we 



use the RAMBO algorithm by Kleiss, Stirling and Ellis [|I^ for this purpose 



An algorithm to insert an additional "soft" particle into a given "hard" configura- 
tion. This algorithm is originally due to Kosower and has been implemented into 



the NLO 4-jet program MERCUTIO g^. 



An algorithm to generate a sequence of events according to a multi-dimensional 



probability density. We use the Metropolis algorithm for that ||17 |. 
In the next two sections we will briefiy review these tools. 



3 Review of the structure of infrared divergences 

In this section we briefiy review the dipole formalism and the factorization properties of 
QCD amplitudes in the soft and collinear limit. We follow the notation of Catani and 
Seymour [ pT^ . 

The NLO cross section is given by 

a^^o = I da"" + I da^ + j da"". (4) 

n n n+1 

da^ represents the LO-contribution, da^ the virtual corrections and da^ the real emission 
part. Taken separately, both da^ and da^ are divergent, only their sum is finite. Within 
the subtraction method, one adds and subtracts a suitable chosen term. The NLO cross 
section is therefore rewritten as 




= j da'' + I \ da^ + I da''\ + I [da"" - da"") . (5) 

n 

The approximation da^ has to fulfill the following requirements: 



j + j [da^ - da^) . 

/ n+1 
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da"^ must be a proper approximation of da^ such as to have the same pointwise 
singular behaviour (in D dimensions) as da^ itself. Thus, da^ acts as a local 
counterterm for da^ and one can safely perform the limit e ^ 0. This defines a 
cross-section contribution 



in^l} - J (d^X^O-d^leJ- (6) 
n+1 

Analytic integrability (in D dimensions) over the one-parton subspace leading to 
soft and coUinear divergences. This gives the contribution 

<f = [ da^+ [ daA . (7) 




The final structure of an NLO calculation is then 

„NLO R I ^NLO I „NLO 



e=0 



All contributions on the r.h.s of eq.(|^) are now finite. However, it should be noted that 

^NLO 
{n+l} 



^ difference of two terms and therefore not guaranteed to be positive-definite. 



Similar, o'^J^'^ contains the interference term between the one-loop amplitude and the 
Born amplitude and is therefore also not guaranteed to be positive-definite. 

The (n + l) matrix element is approximated by a sum of dipole terms 

pairs i,j k^i,j 

= Y Y~T~ — (1, fa'). ...\ ^^ J'^ %fc|l, (ij), fc, •••), 

■ ■ , I - ■ ^Pi ' Pi ii 

pairs t,j k=i't,j * J y 

(9) 

where the emitter parton is denoted by ij and the spectator by k. Here Tj denotes the 
colour charge operator for parton i. The colour charge operators Tj for a quark, gluon 
and antiquark in the final state are 



quark 
gluon 
antiquark 



{...q,...m\...q,...), 
{...g^...\tr%..g\..), 

(10) 



Vij^k is a matrix in the helicity space of the emitter with the correct soft and coUinear 
behaviour. |1, (ij), k, ...) is a vector in colour- and hehcity space. It is also convenient 
to introduce the variable 

Viik = (11) 

2piPj + 2piPk + 2pjPk 

and the ordered set S consisting of all relevant ?/jj^fc-variables: 

S = emitter; j emitted particle; spectator} . (12) 
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We denote the number of elements of S by m. Explicit formulae for the dipole terms 'Dij^k 
have been given by Catani and Seymour |]TB|]. In the CDR-scheme they read [| 



qigj,k 



qiqj,k 



1- Zi{l- Uij^k] 



- (1 + Zi) - eil - Zi 



Pipj 



(ZiP^ - Zjpf){Zip'^ - ZjpJ) 



1 - Zi{l - yij^k) 1 - Zj{l - yij^k) 



+ (1 - e) {zip^ - Zjp^j){zip'( - Zjp]) 

PiPj 



(13) 



where s, s' are the spin indices of the fermion ij and /i, u are the spin indices of the gluon 
ij. The momenta of the emitter and the spectator are related to the original (n + l)-parton 
configuration by 



Pk 



p\ 



1 - yij,k 
Pt + P 



Pk^ 



T, = Pi + Pi - Y——Pk- 



The variables Zi and Zj are given by 



PiPk 



, Zj 1 Zi. 



PjPk + PiPk ' 

We subtract the approximation only for ^ < yj.es- Therefore 
j da"" - da"" 

= j d(pn+l\Mipi,...,Pn+l)fQia\lyij^k>yres)QT+liP^^---^Pn+l) 

+ / d(f)n+l 0-^(^1, ■■■,Pn+l)f (1 - ©(all yij,k > yres)) 0^+*i(Pl' -■^Pn+l) 



(14) 



(15) 



$Z 5Z '^ij'kiPu ■■■,Pn+l)Qiyij,k < 2/res)0r*(Pl' ••• ; Pijj ■■■,Pk, •••) Pn+ly 
pairs i,j k^i,j 



(16) 



Here we split the contribution from da^ into a "resolved" region (first line) and an "unre- 
solved" region (second line). The first line is infrared finite and positive definite. The sum 
of the second and third line is by construction infrared finite, but not necessarily positive 
definite. Both da^ and da^ are integrated over the same {n + 1) parton phase space, but 
it should be noted that da^ is proportional to 6^+i, whereas da^ is proportional to 9™*. 
Here G^"* denotes the jet-defining function for n-partons. 



^Corresponding formulae also exist for the four-dimensional scheme as well as conversion rules between 
the two schemes [|8|. 
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The subtraction term can be integrated over the one-parton phase space to yield the 
term 



I (g) dcr 



B 



j da-^eiyij^k < Vres) = 

-, pairs i,j kj^ij 



Mipole 



T^ij,kQ{yij,k < Vres)- (17) 



The universal factor I still contains colour correlations, but does not depend on the unre- 
solved parton j. The explicit results corresponding to the restricted region of integration 
yij,k < Vres havc bccu givcu by Nagy and Trocsanyi 



dlpdipole^ijfk^i^yijfk ^ Vres) 

1 -e) \2pijpk) 



with 

^qg Vres) 
"^gg Vres) 
"^igi^y Vres) 



271 T( 



rp2 



M 



(n) 
Born 



(18) 





1 






e 




11 








( 


e 



TX 



5-y + 0(e)j, 



9 



2 /I 



3 Ve 
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(19) 



The term I ® da^ lives on the phase space of the n-parton configuration and has the 
appropriate singularity structure to cancel the infrared divergences coming from the one- 
loop amplitude. Therefore 



d(j^ + I O rfa^ 



(20) 



is infrared finite. 



We now introduce a function Qah,c with Qah.c = 1 if Vah.c is the smallest element in 
the set iS, and Oab.c = otherwise. If two invariants are equal (which is for example the 
case for yah,c and yha,c) we consider the invariant yah.c with the smaller value of 



2vhVc 



2paPb + 2paPc + 2pbPc 
to be smaller. Since one of the elements must be the smallest, we have 



ab,c- 



(21) 



(22) 



The sum is over all elements in the set S. We further define a parameter Smin = yminQ"^ 
and slice the region where yas,b is the smallest element in the set into two parts as follows: 
The first region is defined as the region where Sas > Smin and Ssb > Smin- The second 
region is the complement of the first: Sas < Smin or Ssb < Smin- We now rewrite eq.(|l6D 
as follows: 

j da"" -da^ = j dcPn+i J2 ^-^'b (^1 + ^2 + As) (23) 
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with 



Ai = \M{pi,...,Pn+l)fQ{al\yij^k>yres)Q''nMPl,---,Pn+l), 
A2 ©("^as "^min) ('5sf) "^mm) 

■ [\M{pi, .■■,Pn+l)f (1 - 0(all yij,k > yres)) ^n+liPl, ■■■,Pn+l) 

yPijy ■■■jPk, ■■■,Pn+l^ 

pairs i,j kj^i,j 
■^3 [1 0("5as '5mm)0("5s6 "^mm)] 

■ [\M{pi, ...,Pn+l)f (1 - 0(all > y.es)) 0r+l(Pl, --^Pn+l) 

- T^ij,k{Pl^ ■■■,Pn+l)'S){yij,k < yresWnKPl, ■■■ 

pairs i,j k^i,j 



(24) 



Choosing ymin small enough makes the contribution from the As-term negligible. We will 
therefore not consider this term further. From eq. (|23| ) and eq. (^) the structure of our 
algorithm emerges. If all yij^k are larger than yres-, we are in the region corresponding to 
the ^i-term and we generate configurations according to the unsubtracted matrix element 

M.^Bom ■ If one yij^k is smaller than yres, we evaluate the integral of the unresolved region 
corresponding to the y42-term, combine it with the virtual corrections and generate a n- 
parton configuration according to this result. 



4 Review of Monte Carlo tools 

In this section we review some algorithms concerning the generation of phase space as 
well as the Metropolis algorithm for sampling a given distribution. 

4.1 The Metropolis algorithm 

In practical applications one often wants to generate random variables according to some 
probability density P{xi, ...,Xd), which not necesarrily factorizes. In practice one often 
uses the Metropolis algorithm ([|l7l, [T^) for this purpose. Let us call the vector (p = 



{xi,...,Xd) a state of the ensemble, which we want to generate. Within the Metropolis 
algorithm one starts from a state 00, and replaces iteratively an old state by a new one, in 
such a way, that the correct probability density distribution is obtained in the limit of a 
large number of such iterations. The equilibrium distribution is reached, regardless of the 
state one started with. Once the equilibrium distribution is reached, repeated application 
of the algorithm keeps one in the same ensemble. In short, the desired distribution is 
the unique fix point of the algorithm. Two important conditions have to be met for the 
Metropolis algorithm to work: Ergodicity and detailed balance. Detailed balance states 
that the transition probabilities W{(j)i 02 ) and W{(j)2 — ^ <Pi) obey 

P(0l)l^(01 02) = P(02)W^(02^0l). (25) 

Ergodicity requires that each state can be reached from any other state within a finite 
number of steps. Given a state 0i, one iteration of the Metropolis algoritm consists of 
the following steps: 
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1. Generate (randomly) a new candidate 0'. 



2. Calculate A5 = -ln(P(0')/P(0i)). 

3. If AS < set the new state 02 = 4>'- 

4. If AS > accept the new candidate only with probability P{(f)')/ P{(j)), otherwise 
retain the old state 02 = 0i- 

5. Do the next iteration. 

Step 3 and 4 can be summarized that the probability of accepting the candidate 0' is given 
by W{(j)i 0') = min(l, e"^*^). It can be verified that this transition probabilty satisfies 
detailed balance. The way how a new candidate 0' suggested is arbitrary, restricted only 
by the condition that one has to be able to reach each state within a finite number of steps 
as well as the condition that the probability v4(0i cj)') (the probability of suggesting 0' 
given we are in state 0i ) has to be equal to the probability A(0' 0i). 



4.2 Rambo 



The RAMBO-algorithm |TH] maps a hypercube [0, 1]^" of random numbers into n physical 
four-momenta with center-of-mass energy y/P^. Massless four-vectors can be generated 
with uniform weight. 



Let P = [P, 0, 0, 0) be a time-like four-vector. The phase space volume for a system 
of n massless particles with center-of-mass energy y/P^ is 

= /n|^^(^°)^(p')(2^)'^'(^-& (26) 
To derive the RAMBO-algorithm one starts instead from the quantity 

j^xf{x)dx\ . (27) 

The quantity i?„ can be interpreted as describing a system of n massless four-momenta 
that are not constrained by momentum conservation but occur with some weight function 
/ which keeps the total volume finite. The four- vectors are then related to the physical 
four-momenta by the following Lorentz and scaling transformations: 

Pi = X (^75° + b- q?j , Pi = x(qi + bq^ +a(b-qi^b 



where 



n 



i=l 
10 



/ ~ 1 VP2 
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Denote this transformation and its inverse as follows 



b 

By a change of variables one can show 



p'; = xH.^{q,), gf = ii/^fe). (30) 



(31) 



With the choice f{x) = e ^ the integrals over b and x may be performed and one obtains 

Rn = ^n-Sn (32) 

with 

, r (I) r(n- i)r(2n) ^ ^ 

r(n+ 2) 

This gives a Monte Carlo algorithm which generates massless four-momenta pf according 
to the phase-space measure (|26|) . The algorithm consists of two steps: 



1. Generate independently n massless four-momenta gf with isotropic angular distri- 
bution and energies g° distributed according to the density q^e~^^dq^. Using An 
random numbers Ui uniformly distributed in [0, 1] this is done as follows: 



Ci = 2ui^ - 1, ipi = 27rMi2, = - ln(Mi3M 



«4 J ; 



= C[i\/'^ -c^i COS ip,, qf = -c- siny^i, q^ = q°Ci. (34) 

2. The four- vectors are then transformed into the four- vectors p^, using the trans- 
formation (123). 

Each event has the uniform weight 

- - P-)""(i)""rRi^- P^' 

In summary, the RAMBO algorithm provides a mapping from a 4n-dimensional hypercube 
into n physical four- momenta, each event being generated with weight w. By picking 
points in the hypercube at random, we sweep out the complete phase space. Therefore 
ergodicity is satisfied and we can combine the RAMBO algorithm with the Metropolis 
algorithm. 
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4.3 Generating configurations close to soft or collinear regions 

In this section we review an algorithm on how to obtain a (ra + l)-parton configuration 
from a given n-parton configuration. An additional four-momentum ps is inserted be- 
tween the legs Pa and pb- For efficiency we generate the events such that Sas, Ssb > Smin 
with Smin = VminQ^- Wc assume that the complementary region Sas < Smin or Ssb < Smin 
gives only a contribution of order Umin- We first relate a given {n + l)-configuration to 
a n-parton configuration and invert then this relation to obtain an algorithm to generate 



soft or collinear insertions [IG 



Let k'^, kg and k'j^ be the momenta of the (n + l)-parton configuration such that Sas = 
{k'g^ + ksY, Ssb = (^b + ^s)^ and Sab = (^a + ^^ + ^b)^- We want to relate this (n + 1) particle 
configuration to a nearby "hard" n-particle configuration with {ka + kbY = {k'^ + kg + klY, 
where ka and kb are the corresponding "hard" momenta. Using the factorization of the 
phase space, we have 

d$„+i = d^^^,^d^s{K,k',,ks,k',). (36) 
The three-particle phase space is given by 

d^3{K,k'^,ks,k[) = J" — dSasdSsbdQ'^d(j)s 

62{27r)''Sab 

= ■ — dsasdssbd(f)sd^2{K,ka,kb) (37) 

4(^Z7rj Sab 

and therefore 

,^ ,^ dsasdssbdcps 

= 4(2vr)3.., • ^''^ 

The region of integration for Sas and Sgb is Sas > Smin, Ssb > Smin (since we want to restrict 
the integration to the region where the invariants are larger than Smin) and Sas + Ssb < Sab 
(Dalitz plot for massless particles). It is desirable to absorb poles in Sas and Sgb into the 
measure. A naive numerical integration of these poles without any remapping results in 
a poor accuracy. This is done by changing the variables according to 

/ •^mm \ / •-'mm \ /on\ 

Sas = Sab , Ssb = Sab , (39) 

\ Sab J V ^ab J 

where < mi,M2 < 1- Note that ui,U2 > enforces Sas,Ssb > Smin- Therefore this 
transformation of variables may only be applied to invariants Sij where the region < 
Sij < Smin is cut out. The phase space measure becomes 

d^ri+l = d^n ,,^ SasSsb f Snun\ q^^^^ ^ ^ Sab)dUidU2d(f)s- (40) 

4(27r)'^ Sab V ^ab ) 

This give the following algorithm for generating a (n + l)-parton configuration: 
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1. Take a "hard" n-parton configuration and pick out two momenta ka and ki,. Use 
three uniformly distributed random number Vi,i'2,V3 and set 



Sab ^ 


S-min 


Sab 


Sob ^ 


Smin 


Sab 


27rvs 





V2 



(41) 

2. If {sas + Ssb) > Sab, rejcct the event. 

3. If not, solve for k'^, k'^ and kg- If Sas < Sgb we want to have k^ kb as Sas 0. 
Define 

J-, Sab Ssb j-i Sas ^" Ssb j-. Sab Sas //loN 

^a = o ; ^ -^s = o / > -^6 = c) , 1 i^^J 

^ySab ^ySab Sab 



9ab — arccos ( 1 — — — — — ) , Ogb — arccos ( 1 I • (43) 

V '^EaEb 7' V '^EsEbJ ^ ' 

It is convenient to work in a coordinate system which is obtained by a Lorentz 
transformation to the center of mass of ka + kb and a rotation such that k'^ is along 
the positive z-axis. In that coordinate system 

= Ea(l,sin6'a6cos(0s + 7r),sin6'a6sin(0s + 7r),cos6'ab), 
Ps = E5(l,sin6's6cos0s,sin6'sbsin03,cos6'sb), 

= ^6(1, 0,0,1). (44) 

The momenta p^, "Ps and p'j, are related to the momenta A;„, kg and A;^ by a sequence 
of Lorentz transformations back to the original frame 

K = hoostKy{(t>)KzWa (45) 

and analogous for the other two momenta. The explicit formulae for the Lorentz 
transformations are obtained as follows : 



Denote by 7^ = ^/(k^^^+~kbf and by pb the coordinates of the hard momentum 
kb in the center of mass system of ka + kb- Pb is given by 

with Q = ka + kb- The angles are then given by 



= arccos Tl-^l 



<t> = arctan(^^). (47) 
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The explicit form of the rotations is 



The boost k' = IVboostfl is given by 



/ 


1 





\ 







cos^ 


sin^ 




U 


U 


i U 


\ 


u 


— sin 6 


U COS U J 


/ 


1 





\ 







cos 


— sin 







sin (p 


COS0 


V 








l) 



k' 



q-Q 



K(QO + K) K 



with Q = ka + h and K = -\/(^a+"^bF- 

4. If Sas > Ssb, exchange a and b in the formulae above 

5. The "soft" event has then the weight 

1 SasSsfe 



TT 



Wn+l 



2(27r)3 Sab 



In' 



Sab 



Wr, 



(48) 



(49) 



(50) 



where is the weight of the original "hard" event. 
Again by picking vi, V2 and ^3 randomly, we sample the complete phase space except the 

region Sas Smin Sgb ^ Sjnin- 



5 Generating unweighted events according to leading- 
order matrix elements 

In this section we give an algorithm how to generate a sequence of unweighted events 
according to leading-order matrix elements. This section serves as a warm-up exercise 
to the next section, where we extend the algorithm to next-to-leading order matrix el- 
ements. Since for leading-order calcultions the amplitude squared is always positive, 
negative weights do not occur at this stage. However wc would like to keep the discussion 
quite general and process- independent. Wc therefore do not make any assumptions on 
the form of the matrix element squared. In particular we do not assume that the am- 
plitude squared, viewed as a multi-dimensional probability density factorizes nor do we 
assume that we have at our disposal a simple function which is everywhere larger than 
the amplitude squared. Therefore the acceptance-rejection method is not at our disposal 
and we use the Metropolis algorithm to generate the desired distribution. 

Observables in high-energy collider experiments are often of the following form 

(O) = j d^n{Pa+Pb:Pl:-:Pn) J{0,Pu...,Pn): (51) 
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where pa and pb are the momenta of the incoming particles, the outgoing particles are 
denoted by the labels l,...,n. The Lorentz-invariant phase space is denoted by d^n, 
l/{8K{s)) is a kinematical factor with s = {pa + Pb)"^ which includes the averaging over 
initial spins (we assumed two spin states for each initial particle), \A4f is the relevant 
matrix element squared and J{0,pi, ...,Pn) is a function which defines the observable and 
includes all experimental cuts. The situation described above corresponds to electron- 
positron annihilation. If hadrons appear in the initial state there are slight modifications. 
Using the RAMBO algorithm eq. (|5TD may be rewritten as 

(O) = I rf^"tx^o(^) '^g^g|y' J(0,pi,...,p„), (52) 

where the weight wo{u) is given by eq. (|35|) and the final-state four-momenta are obtained 
from the Mj's through the RAMBO algorithm. Suppose now that we have at our disposal a 
jet-defining function Qy^^^iPi, ■■■,Pn) such that the support for J{0,pi, ...,p„.) is contained 
in the support for ©^^^^(Pi, ■■■,Pn)- The function ©^^^^(Pi, ■■■,Pn) returns 1 if the n partons 
are resolved at the scale i/res and otherwise. We may then replace J(0,pi, ...,Pn) by 

eiJpi,...,Pn)J{0,p,,-.,Pn) (53) 

and generate random numbers uj distributed according to the probability density function 

\M{p^{u)fQlJp,,...,Pn), (54) 
(Jo 8K[s) ^ 

where cxo is the total LO-cross section, which normalizes the probability density function 
to unity: 

^0 = ^^^ld'-uwo{u)\M{p.{u)fQlJp^,...,p^). (55) 

Note that the quantity in eq. (Q) is positive-definite, as required for a probability density 
function. We can apply the Metropolis algorithm to generate the desired distribution f\. 

1. To suggest a new candidate we randomly pick a point {u[, u'^^ in the 4n-dimen- 
sional hypercube and use the RAMBO algorithm to obtain the corresponding n- 
parton configuration of four- vectors p'^, p^. 

2. The probability density P{u') for this candidate is given by eq. (|5^. 

3. If P{u') = (which can occur due the presence of the jet-defining function 
^yres(p'i^ ■■■yP'n)) reject the candidate and go back to step 1. 

4. Calculate AS* = — ln(P(n')/P(n)), where P{u) is the probability density of the 
previous state. 

5. If AS < accept the new candidate and return the n-parton configuration p[, 
P'n- 

^The Metropolis algorithm has also been considered by Kharraziha and S. Moretti for the gen- 
eration of phase space. These authors use a different method to suggest a new candidate and do not 
combine the Metropolis algorithm with a "standard" algorithm to generate the phase space. 
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6. If AS > accept the new candidate only with probabihty P{u')/ P{u), otherwise 
retain the old state and return the old n-parton configuration pi,...,Pn. 

7. Do the next iteration. 

A technical comment on step 3 is in order: Step 3 immediately rejects candidates which do 
not pass the cuts O^,,^^ • This does not violate detailed balance due to the specific way new 
candidates are suggested in step 1: A new state (m'^, ^4^) is suggested independently 
of the present state (mi, ...,M4„). If one uses a different procedure in step 1, for example 
by changing the k-th coordinate Uk ^ Uk + S hj a small, but random value S, candidates 
suggested from states close to the bounday of B^^^^ are more likely to fall outside the 
boundary then candidates suggested from states inside the region. To ensure detailed 
balance, one has to discard step 3 and treat candidates with zero probability in step 6 as 
candidates which are never accepted. This inserts the old state once again into the event 
record. 



6 Generating unweighted events according to next- 
to-leading-order matrix elements 

In this section we generalize the algorithm from the previous section to next-to-leading 
order. For the real emission part we start from eq. p3|) and eq. (0) and first describe 
a method, how to pick two legs a and b corresponding to the emitter and the spectator. 
We then try to insert a soft particle in between these two. We divide the interval [0, 1] 
into m subintervals (m is the number of elements in the set S) and denote the size of the 
l-th sub interval by Pi. We further introduce functions Qi{x) which return 1 if x is in the 



l-th. subinterval and zero otherwise. We may therefore rewrite eq. (|23|) as 

J da"" -da^ = JdxJ2 J5-^^ii^) j d(pn+iQas,b (A, + A2) . (56) 

Now the value of x decides where the soft leg gets inserted. In more detail, we may 
rearrange the above equation as 

Wo{u) j dx ^ -^0/(2^) j dVi dV2 dv3 WsofSas,b {M + A2) , 

(57) 



da^ - da' 



with 

^ 1 ^as^sb I 2 I ^min\ /ro\ 

= 2(2^— (i:rj- 

We may therefore use 4n -|- 1 -|- 3 random numbers mi,...,M4„, x, f i, f 2 and f 3, all uniformly 
distributed in [0, 1] and generate the phase space as follows: We use the random numbers 
Mi,...,M4„ to generate a hard configuration of n four- momenta using the RAMBO algo- 
rithm. The value of the random number x decides where we try to insert a soft leg. We 
check into which of the m subintervals x falls and pick out the corresponding element of 
the set S. This tells us which parton is going to be the emitter and which parton is going 
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to be the spectator. We now use the random numbers t'l , t'2 and vs and the algorithm for 
soft insertions and generate a (n + l)-parton configuration, where a soft particle has been 
inserted between the emitter and the spectator. This suggests the following Metropolis 
algorithm: 

1. To suggest a new candidate we randomly pick a point {u'l, u'^^, x', v[, V2, v'^) in the 
4n + 1 + 3-dimensional hypcrcube. Using the first An random numbers ui,...,U4n and 
the RAMBO algorithm we obtain the corresponding "hard" n-parton configuration 
of four- vectors (^'f^^rd — ?n)- Ftom this hard configuration and the values 
of X and Vi, V2 and we obtain a related "soft" {n + l)-parton configuration 
^'soft ~ (p'l' P'n^P'n+i) together with the information where the soft particle s has 
been inserted. Therefore the emitter a and the spectator b are known. 

2. Check if yas,b is the smallest value in the set S. If not, reject the candidate and go 
back to step 1. 

3. Check if all yij^k are larger then yres- 



If this is the case, the probabilty density of the new candidate is given by 



P' = 



1 wo{u)w 



soft 



(Tnlo 8K{s) Pi 
and we set a flag f — 1- Here (Tnlo is given by 
1 



cnlo 



SK{s) 



M 



(n) 
Born 



+ 2 Re Mitn Ml\o, + {M^Znim'Zn) 



(n) 



6" 

yres 



+ / d(pn+l 



M 



(n+l) 
Born 



yres 



■ (60) 



pairs i,j k^i,j 

Oy^2 returns 1 if at least n partons are resolved at the scale yres- 
If P' = (which can occur due the presence of the jet-defining function 
,Pn+i)) we reject the candidate and go back to step 1. 



If at least one yij k in the set S is smaller then y^es the probability density is given 

by 

1 I Wq 



GNLO SK{s) V 



M 



(n) 
Born 



+ 2 Re A^S.;A^tU + {MZnimZn 



(n) 



(61) 



where the integral over the unresolved region is given by 



R = ^ M^"! dV2 dV3 WsofSas,b 



M 



(n+l) 
Born 



[l-Q{Ql\yi^,k>yres))% 



n+l 
'res 



pairs i,j k^i,j 



(62) 
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and the unresolved volume is given by 



V 



^ / dVi dV2 dV3 Qas,b (1 - ©(all yij,k > Vres)) 



(63) 



s 



In addition set the flag / = 0. 

If P' < (which can occur if y^es is chosen too small) throw an exception. 
4. Calculate AS* = — ln(P'/P), where P is the probability density of the previous 



5. If AS < accept the new candidate. If the flag / = 1 return the "soft" (r;,+l)-parton 
conflguration p'l, p'n otherwise if / = return the "hard" n-parton conflguration 

q'i,--,qn- 

6. If AS > accept the new candidate only with probability P'/P, otherwise retain 
the old state and return the old parton conflguration. 

7. Do the next iteration. 

If the value of y^es is chosen too small, it can occur that P' in eq. (^) turns out negative. 
However, this is not an artefact of our algorithm, but a failure of perturbation theory: P' 
corresponds to the NLO cross section for observing n jets with (flxed) four- vectors q[, 

resolved at a scale i/res- If Vres is too small the NLO prediction might become negative. 
This is of course unphysical. The algorihtm throws an exception in this case and one 
should rerun the program with a larger value of t/res- 

Given a generated sequence of momentum conflgurations 0i,...,0Ar, where each momen- 
tum conflguration contains either n or (n + 1) four- vectors, the Monte Carlo estimate for 
an infrared safe observable is then given by 



Here J{0, (pi) is a function which deflnes the observable and includes all experimental 
cuts. The approximations involved are: 

• If within an {n + l)-parton conflguration some value of ytj^k is smaller then t/res, this 
conflguration is replaced with an n-parton conflguration, which corresponds to the 
average over the unresolved region combined with the virtual corrections and the 
Born term. 

• Due to the presence of the jet-deflning functions ©^^^^ and ©^^^ conflgurations 
which would classify as (n — l)-parton events according to these functions are cut 
out. One has to ensure that the support for the actual observable O is contained in 
the support for 9" and ©I^"*"^. 

• We neglected the term A^, in eq. (^). This introduces a (negligible) systematic 
error of i/min- This is done for efficiency. 



state. 




(64) 



1=1 
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7 Conclusions 

In this paper we gave an algorithm to generate unweighted events of n or (n + 1) four- 
vectors according to next-to-leading order matrix elements. We considered massless QCD 
in electron-positron annihilation. The extension of our methods to initial-state partons 
or massive partons should be possible. 
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